Tomography refers to the cross-sectional imaging of an object from either transmission, emission or reflection data collected from many different directions. Tomographic imaging deals with reconstructing an image from such data, commonly called projections. From a purely mathematical standpoint, a projection at a given angle is the integral of the image in the direction specified by that angle. Most of the powerful new medical imaging modalities that have been introduced during the last four decades, such as Computed Tomography (CT), Single Photon Emission Computed Tomography (SPECT), Positron Emission Tomography (PET), Magnetic Resonance Imaging (MRI) and 3D ultrasound (US), were the result of the application of tomographic principles.
The present invention relates to general tomographic reconstruction problems and, therefore, can be adapted to address all the aforementioned imaging modalities. However, for convenience and simplicity, the present specification is concerned with resolving PET image reconstruction problems, keeping in mind that the present invention could be used to resolve general tomographic reconstruction problems.
Positron emission tomography is a non-invasive imaging modality which is said to be metabolic or functional since it allows to measure and localize the biodistribution of a radio-element injected inside the body under investigation. The fundamentals of PET imaging is based on the localization of a positron source within the tissues of the subject by the detection of two oppositely directed gamma rays, emitted during the annihilation of the positron. An approximate localization of the site of disintegration along a line is possible by detecting simultaneously the two annihilation photons with opposed detector of the camera. This region where the disintegration could have taken place, known as “line of response (LOR)”, “tube of response (TOR)”, or “coincidence line”, is defined by the volume joining the two detectors that received the annihilation photons. The time of detection of both annihilation photons must be measured and compared to validate that the two photons where emitted at the same time and thus, have a high probability of coming from the same disintegration. A PET scan consists of measuring and comparing the time of occurrence of every detected photon in the camera to record the coincident events that occur between all possible detector pairs which represent all the TORs of the camera. Using those measurements, an image of the biodistribution of the tracer injected in the subject could be obtained with image reconstruction methods.
A PET imaging system is typically composed of many 511 keV gamma ray detectors disposed on at least one ring or on a stack of many rings of detectors. The system also include a signal processing unit used to extract information relevant to the PET measurements (e.g. the time, the energy and the position of the detected photon) and a coincident sorter unit used to verify and count coincident events between all pairs of detectors in the camera. The imaging system can be used in two-dimensional acquisition mode where coincidence measurements are performed only between detectors in the same ring. Some apparatus also allow three-dimensional acquisition mode where coincidence measurements can be performed between detectors at different ring positions.
The geometries and the disposition of all sensing elements in a system is an important aspect when addressing the problem of reconstructing an object from a measure of its projections. A first solution to this problem date back to the paper by Radon in 1917 which proposed the inverse-Radon transform. Those works lead to the development of analytical reconstruction algorithms which aim at solving the inverse-Radon transform by considering every tubes of response (TORs) of the camera as an infinite thin line joining the two coincident detectors. One of the most common algorithm of this class is the filtered backprojection (FBP) algorithm, which makes use of the projection-slice theorem. Some drawbacks of the method is the presence of artifacts and the misplacements of some radioactive objects in the reconstructed image due to the over-simplification of the TOR functions by an infinite thin line.
A solution to this problem consists of modeling more closely the function representing the detection sensibility along the TORs which depends on the geometries and on the position of the two coincident detectors. Such a function is commonly called a coincidence aperture function (CAF) and could be obtained analytically or from a Monte Carlo simulation. Using the CAF of every TOR of the camera, a probability matrix relating the radioactive density of every pixel in the image to the measurements made on every TOR of the camera could be obtained. The problem could then be represented by the following relation:y=Af  (1)where each coefficient ai,j of the probability matrix A represents the probability that an event produced in the jth pixel of the image vector f has been detected by the ith detector pair of the measurement vector y. The resolution of this problem could be addressed by two broad class of image reconstruction methods: “direct methods” which consist of inverting the probability matrix to obtain the image directly from a vector-matrix multiplication between the inverted matrix and the camera measurements and “iterative methods” which consist of making some successive estimates of the density distribution of the image in respect to the probability matrix and to the measurements until the image converge to a solution that meet given criteria.
Direct and iterative methods both have some advantages and drawbacks. Irrespectively to the method used, the main goal of reconstructing an image from its projections is to obtained an image that is as close as possible to the true density distribution of the object being imaged. In that respect, characteristics like spatial resolution, image contrast, signal-to-noise (S/N) ratio are aspects which can improve the detection accuracy but these should not be obtained at the expense of having some artifacts in the image which could lead to false detections. The computation time of the algorithm can also play a role in the selection of the most appropriate algorithm to use for a given scanner and in a given context of utilization like in clinical applications. A review of prior methods will be made in the light of those considerations.
Direct Methods:
A direct image reconstruction method based on pseudo-inversion of matrices may be found in [Llacer, Tomographic image reconstruction by eigenvector decomposition: Its limitations and areas of application]. The idea was to multiply both side of equation 1 by the transpose of the probability matrix which result in:ATy=ATAf  (2)and to obtain the image by inverting the matrix ATA which can be expressed as:f=(ATA)−1ATy  (3)
The matrix to inverse is now ATA and the pseudo-inversion operation is facilitated by the fact that ATA is a symmetric matrix. In this form, the matrix ATA can be viewed as a two-dimensional blurring matrix of the PET imaging system since it is obtained from the product of a backprojector AT and a projector A of the imaging system. One of the drawback of the method is the high condition number (CN) (largest divided by smallest eigenvalue) of the blurring matrix for tomographic reconstruction problems. The high CN leads to imprecision in the computation of eigenvectors corresponding to the smallest eigenvalues. Another drawback is the computational burden of the pseudo-inversion operation. Moreover, the measurement must be backprojected by AT before being multiplied by the pseudo-inverse of ATA which, however, does not add much computation compared to the pseudo-inversion operation.
A method to accelerate the blurring matrix ATA pseudo-inversion was proposed by Barber [Barber, Image reconstruction, “European patent specification”, 1992]. The method consists of taking advantage of the block circulant structure emanating from the blurring matrix ATA when the projector A and the backprojector AT are strip functions defined in the continuous domain. The pseudo-inversion of the block circulant blurring matrix can then be performed independently on smaller sub-matrices by using the properties of the Fourier transform. The discretization of the image into square pixel elements is performed during the last step of backprojecting the result of the multiplication between the inverted blurring matrix and the measurements which lead to the following relation:f=AT(AAT)−1y  (4)
One can notice that equation 4 is equivalent to equation 3, but with a different operation ordering. Barber has also pointed out that the backprojector can be some other function than the transpose of the projector.
The coefficients of the block circulant ATA matrix can be viewed as a natural pixel decomposition of a two-dimensional image. Buonocore [Buonocore, A natural pixel decomposition for two-dimensional image reconstruction, 1980] was the first to propose the term “natural pixel” to define pixels arising naturally out of the geometry of the X-ray beam paths used to measure the projections. One advantage of such an image decomposition is that it preserved the symmetries between the TORs of a camera leading to a block circulant structure of the ATA matrix. The projector A and the backprojector AT used in this formulation of the problem are strip functions defined in the continuous domain. The discretization of the forward projection matrix A into square pixels in a Cartesian grid representation, like the projector used in [Llacer, Tomographic image reconstruction by eigenvector decomposition: Its limitations and areas of application], would broke the symmetries between the TOR functions and will therefore not lead to a block circulant blurring matrix. This observation was pointed out by Baker who also proposed using natural pixel to facilitate the inversion problem in tomographic image reconstruction [Baker, “Generalized approach to inverse problems in tomography: Image reconstruction for spatially variant systems using natural pixels”, 1992]. In this work, Baker proposed to accelerate pseudo-inversion of the blurring matrix through the used of the singular value decomposition (SVD) on a Fourier transformed version of the block circulant blurring matrix. To obtain the final image, the result of the multiplication between the pseudo-inverse blurring matrix and the measurement vector is backprojected.
It was shown by Shim [Shim, “SVD Pseudoinversion image reconstruction, 1981] that the pseudo-inverse of matrices by singular value decomposition techniques can be performed directly on the system probability matrix A, as stated in equation 1. One advantage of performing the SVD directly on the system matrix A is that the CN for such problems will be lower then one obtained for the blurring matrix ATA. When a lower CN is encountered, the SVD algorithm will be less affected by computer precision limitations which will lead to better estimate of smaller singular values and of the corresponding singular vectors. Another advantage of the method is that the reconstructed image can be obtained directly from the multiplication of the inverted matrix with the measurement vector which avoids the used of a backprojection matrix. In other words, the quality of the reconstructed image is not influenced by the modeling of the backprojection matrix. This is a significative improvement since it is easier and more precise to include all the physical aspects of the imaging process in the projection matrix than in the backprojection matrix. Vandenberghe [Vandenberghe, “Reconstruction of 2D PET data with Monte Carlo generated natural pixels] has shown that it is possible to include a more realistic modeling of the imaging process inside the natural pixels decomposition emanating from the blurring matrix ATA by using a Monte Carlo simulation. However, Vandenberghe fails to include the modeling of the system in the backprojection matrix.
A drawback of performing the SVD directly on the probability matrix A is the computational burden of decomposing such a large matrix for ill-posed problems. Selivanov [Selivanov, “Fast PET image reconstruction based on SVD decomposition of the system matrix] proposed to perform the SVD step only once for a given PET scanner configuration. The resulting image reconstruction algorithm is really fast since it required only a matrix-vector multiplication between the inverted probability matrix and the measurement vector. However, due to limitations in computation power, available memory size and floating point precision of current computer, the SVD inversion can only be performed for modest size imaging problem.
Iterative Methods:
Iterative image reconstruction methods consist of estimating the density distribution of the image in respect with the system probability matrix and the measurements made on all TORs of the imaging system. Two wide class of methods are based on this approach: Algebraic reconstruction methods and statistical reconstruction methods. Algebraic reconstruction methods consist of resolving directly the system of equation by some iterative process. Instead of solving directly the system of equations, statistical image reconstruction methods can include some probability distribution in the iterative process. Those methods thus have the advantage of taking into account the stochastic nature of the detection process. This could lead in a better estimate of the object density, especially in the case of low projection statistics.
It has been shown that iterative reconstruction methods, based on a probability matrix of the system, could lead to images with higher spatial resolution and with a more uniform signal-to-noise (SNR) distribution than analytical methods like FBP. Nevertheless, a drawback of iterative methods is the computational burden coming from the many iterations required before the image reach final convergence. Moreover, the number of computation required at each iteration step become more important as the size of the probability matrix A increased (equation 1). The size of the probability matrix depends on the number of TORs measured by the camera and on the number of pixels in the reconstructed image and also on the accuracy of the acquisition model use to derive the matrix. For an imaging system composed of many detector rings and operating in full three-dimensional (3D) acquisition mode, the events detected in millions of TORs can be recorded leading to a very huge system matrix and therefore required a lots of computations when ‘true’ 3D iterative image reconstruction is performed. By ‘true’ 3D iterative image reconstruction methods, we mean a method that use all relevant information collected by the camera to estimate every pixels of the image.
In contrast, some ‘approximate’ 3D iterative image reconstruction methods use data rebinning of 2D reconstructed image to estimate the 3D image distribution. In such algorithm, the 3D projection data is first sorts into smaller 2D data set containing the TOR measurements for each transaxial slice to be reconstructed. All the different slices are then reconstructed independently using a 2D iterative image reconstruction algorithm based on a smaller probability matrix which relates only the measurements and pixels belonging to the same 2D slice. All the 2D slices are then rebinned together to estimate the 3D image distribution. By decomposing the 3D image reconstruction problems into smaller independent 2D image reconstruction problems, a significant reduction of the computation time is achieved. The use of a smaller 2D probability matrix also avoids the cumbersome handling of a huge 3D probability matrix. Nevertheless, this oversimplification of the initial problems can lead to mispositioning and/or to non-optimal estimate of some source activity in the 3D image.
To accelerate ‘true’ 3D iterative image reconstruction methods and make them fast enough for day-to-day use in clinical applications, one can address one or more of the following problems: 1) the number of iterations required by the algorithm to converge, 2) the number of operations required in each iteration loop and 3) the size of the system matrix. The acceleration of 3D iterative methods can also be obtained by sharing the task of computation over many parallel processors and/or by using a more powerful computation unit.
In the literature, most effort in accelerating iterative image reconstruction methods have been concentrated in addressing the first aspect which consists of reducing the number of iterations required to converge to the optimal image. Probably the most known example of such an acceleration technique is the Ordered Subsets Expectation Minimization (OSEM) algorithm [Hudson, “Accelerated image reconstruction using ordered subsets of projection data”, 1994] which leads to acceleration of the convergence of the well known Maximum Likelihood Expectation Minimization (MLEM) algorithm that was first proposed by Shepp and Vardi [Shepp, “Maximum likelihood reconstruction for emission tomography”, 1982]. The OSEM algorithm consists in dividing the projection measurement into small subsets (or blocks) and updating the image estimate using only the data of one subset at a time. One complete iteration loop is completed when all the subsets have been used in the update of the image estimate. This strategy provides an order-of-magnitude acceleration over the MLEM algorithm which is proportional to the number of subsets used. Many other solutions to accelerate the convergence of iterative algorithms are proposed in the literature but will not be presented here since they have less to do with the present invention.
Another strategy for accelerating the speed of iterative algorithms consists of minimizing the number of operations involved in each iteration loop. For most iterative image reconstruction algorithms the main computational burden comes from the matrix-vector operations involved in the forward and back projection steps that are performed at each iteration loop. The number of operations in the forward and back projection steps depends on the size and on the sparsity of the system matrix. In PET, the matrix coefficients are non-null only for pixels (or voxels) which have a non-null probability of emitting a disintegration that is detected by a given TOR. Using a simplified geometrical model of the acquisition process, many matrix coefficients are null and do not have to be stored. By storing the system matrix in a sparse format, both the memory requirements and the computational burden are reduced by using only the non-null values while performing computations in the forward and back projection steps. However, for some complex 3D image reconstruction problems, the sparse probability matrix could still be very huge.
The size of the probability matrix is really becoming a problem when it is too large to fit in the random access memory (RAM) of a given processor or dedicated hardware computational unit. A solution consists in storing the precomputed system matrix in a larger memory location, usually a hard disk, and to access only sub-part of the precomputed system matrix during the forward and back projection steps. The overhead of reading large amount of information on a hard disk could lead to very long delay for the image reconstruction procedure. Another solution consists in recomputing on-the-flag all the system matrix coefficients at each iteration loop to avoid storing them. Some simplifications of the acquisition model should however be used while computing the system matrix coefficients to keep the computation time within acceptable delay. Those simplifications in the system matrix will affect the quality of the reconstructed image.
A strategy for reducing the memory requirements for storing the system matrix consists of taking advantage of the symmetries between the TORs of a camera to store only non-redundant part of the system matrix. In-plane symmetries come from the angular repetitions of blocks of detectors within the ring of an imaging system. In-plane angular symmetries can also be obtained by performing successive measurements taken at different angle position using for example a rotating detector gantry. For camera operating in 3D acquisition mode, one can also take advantage of the axial symmetries present between the different detector rings. Some mirror symmetries could also be retrieved for some camera configurations. Most image reconstruction methods presented in the literature are based on a Cartesian image representation which have the consequence of broking most of the system symmetries in the system matrix. A solution consists of using a basis function defined according to a polar or cylindrical coordinate grid that is specifically designed to preserve all symmetries of a given camera configuration in the system matrix. The idea of using a polar image was first proposed by Kearfott [Kearfott, K. J., “Comment: Practical considerations;”, Journal of the American Statistical Association, March 1985, pp. 26-28]. This solution was proposed to overcome the memory limitation of computer of that time, according to the size of probability matrix used in two-dimensional iterative image reconstruction techniques. One drawback of the polar image configuration used by Kearfott is that a fixed number of pixels are used at every radius position which has the consequence of dividing the image into very thin pixels at the center and into wide pixels at the border of the field of view. The important size disparities between the innermost and the outermost pixels of the image could result in over resolution in the center of the image and in degradation of the imaging system spatial resolution farther from the center. Moreover, Kaufman argued that, for statistical reasons, pixels with similar area should be used by iterative algorithms.
To overcome the problem of size disparities between the polar pixels, Kaufman [Kaufman, Implementing and accelerating the EM algorithm for positron emission tomography] proposed to use a polar image with different number of pixels at each radius position and with variable distances between successive radius positions in such a manner that pixels of similar area are obtained. In this paper, Kaufman showed that important memory saving can be obtained by using a polar image in replacement of a Cartesian image. However, the proposed polar image reconstruction algorithm did not result in significant acceleration of the forward and back projection steps.
Hebert [Hebert, Fast MLE for SPECT using an intermediate polar representation and a stopping criterion] also proposed the use of a polar image representation. In this work, it was shown that the polar image based system matrix can be reordered into a matrix having a block circulant structure. The block circulant system matrix can then be converted in the Fourier domain in order to reduce the number of operations involved in the forward and back projection steps of iterative image reconstruction methods. One drawback of converting the probability matrix in the Fourier domain comes from the fact that a lot of null values in the matrix become non-null during the Fourier transform which increase the matrix size. Another drawback of the method proposed by Hebert, is that the polar image used has a constant number of pixels at every radius position which results in size disparities between innermost and outermost pixels. To avoid the potential problem of using pixels with size disparities in iterative algorithm, Hebert converts the polar image representation into a Cartesian image representation before performing the image update. However, this conversion from a polar to a Cartesian image representation and then from a Cartesian to a polar image representation can result in a loss of spatial resolution due to divergences between the pixels size and position of the polar and Cartesian images. Moreover, since a dense system matrix in the Fourier domain is used in the computation, the gain of speed compared to the traditional method based on a sparse system matrix in the spatial domain rapidly vanishes as the number of symmetries in the camera become small compared to the number of detectors within a ring.